自己相関行列の逆行列の観察¶
なんといってもまずは観察から。
In [1]:
import math
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
NUM_SAMPLES = 4096
# データ生成
sindata = np.zeros(NUM_SAMPLES)
gaussnoise = np.random.normal(0, 1, NUM_SAMPLES)
laplacenoise = np.random.laplace(0, 1, NUM_SAMPLES)
for i in range(NUM_SAMPLES):
sindata[i] = math.sin(2 * math.pi * 440 * i / 48000.0)
In [2]:
# 自己相関行列を計算
def calc_auto_corration_matrix(data, dim):
S = np.zeros((dim, dim))
for smpl in range(dim, data.size, 1):
xvec = data[smpl - dim:smpl].reshape((dim, 1))
S += np.dot(xvec, xvec.T)
return S / (data.size - dim)
# 分散1, 平均0化
def normalize(data):
return (data - np.mean(data)) / np.std(data)
# 自己相関行列とその逆の表示
def plot_gt_auto_correlation(data, dim):
S = calc_auto_corration_matrix(data, dim)
_, axs = plt.subplots(ncols=2, figsize=(10, 4))
sns.heatmap(S, ax=axs[0])
axs[0].set_title('Auto Correlation (Grand Trues)')
sns.heatmap(np.linalg.inv(S), ax=axs[1])
axs[1].set_title('Inverse Auto Correlation (Grand Trues)')
In [3]:
# ノイズの比をいじって調べる
plot_gt_auto_correlation(normalize(gaussnoise), 16)
plot_gt_auto_correlation(normalize(sindata + 1.0 * gaussnoise), 16)
plot_gt_auto_correlation(normalize(sindata + 0.75 * gaussnoise), 16)
plot_gt_auto_correlation(normalize(sindata + 0.5 * gaussnoise), 16)
plot_gt_auto_correlation(normalize(sindata + 0.25 * gaussnoise), 16)
plot_gt_auto_correlation(normalize(sindata + 0.0 * gaussnoise), 16)
In [4]:
# 信号に1次の相関を付与
def make_1st_correlation(data, corr):
corr_data = data.copy()
for smpl in range(1, data.size, 1):
corr_data[smpl] += corr * corr_data[smpl - 1]
return corr_data
corr_gaussnoise = make_1st_correlation(gaussnoise, 0.95)
plot_gt_auto_correlation(normalize(corr_gaussnoise), 16)
plot_gt_auto_correlation(normalize(sindata + 1.0 * corr_gaussnoise), 16)
plot_gt_auto_correlation(normalize(sindata + 0.75 * corr_gaussnoise), 16)
plot_gt_auto_correlation(normalize(sindata + 0.5 * corr_gaussnoise), 16)
plot_gt_auto_correlation(normalize(sindata + 0.25 * corr_gaussnoise), 16)
plot_gt_auto_correlation(normalize(sindata + 0.0 * corr_gaussnoise), 16)
In [5]:
for ratio in np.linspace(0.0, 1.0, 10):
mixed = ratio * normalize(sindata) + (1.0 - ratio) * normalize(gaussnoise)
mixed = normalize(mixed)
plot_gt_auto_correlation(mixed, 16)
for ratio in np.linspace(0.0, 1.0, 10):
mixed = ratio * normalize(sindata) + (1.0 - ratio) * normalize(corr_gaussnoise)
mixed = normalize(mixed)
plot_gt_auto_correlation(mixed, 16)
In [6]:
# 自己相関行列の逆行列の比較: (data+noise, noiseのみ)の比較
def plot_compare_inverse_auto_correlation(mixed, noise, dim):
Smix = calc_auto_corration_matrix(mixed, dim)
Snoise = calc_auto_corration_matrix(noise, dim)
Smix_inv = np.linalg.inv(Smix)
Snoise_inv = np.linalg.inv(Snoise)
_, axs = plt.subplots(ncols=3, figsize=(14, 4))
sns.heatmap(Smix_inv, ax=axs[0])
axs[0].set_title('Inverse Auto Correlation for Data+Noise')
sns.heatmap(Snoise_inv, ax=axs[1])
axs[1].set_title('Inverse Auto Correlation for Noise')
sns.heatmap(np.abs(Smix_inv - Snoise_inv), ax=axs[2])
axs[2].set_title('Absolute Error')
for ratio in np.linspace(0.0, 1.0, 10):
mixed = ratio * normalize(sindata) + (1.0 - ratio) * normalize(gaussnoise)
mixed = normalize(mixed)
plot_compare_inverse_auto_correlation(mixed, normalize(gaussnoise), 16)
for ratio in np.linspace(0.0, 1.0, 10):
mixed = ratio * normalize(sindata) + (1.0 - ratio) * normalize(corr_gaussnoise)
mixed = normalize(mixed)
plot_compare_inverse_auto_correlation(mixed, normalize(corr_gaussnoise), 16)
In [7]:
# 自己相関行列の逆行列の比較: (data+noise, noiseのみ)の比較 L2ノルム正規化して比較
def plot_norm_compare_inverse_auto_correlation(mixed, noise, dim):
Smix = calc_auto_corration_matrix(mixed, dim)
Snoise = calc_auto_corration_matrix(noise, dim)
Smix_inv = np.linalg.inv(Smix)
Smix_inv /= np.linalg.norm(Smix_inv)
Snoise_inv = np.linalg.inv(Snoise)
Snoise_inv /= np.linalg.norm(Snoise_inv)
_, axs = plt.subplots(ncols=3, figsize=(14, 4))
sns.heatmap(Smix_inv, ax=axs[0])
axs[0].set_title('Inverse Auto Correlation for Data+Noise')
sns.heatmap(Snoise_inv, ax=axs[1])
axs[1].set_title('Inverse Auto Correlation for Noise')
sns.heatmap(20 * np.log10(np.abs(Smix_inv - Snoise_inv)), ax=axs[2])
axs[2].set_title('Absolute Error [dB]')
for ratio in np.linspace(0.0, 1.0, 10):
mixed = ratio * normalize(sindata) + (1.0 - ratio) * normalize(gaussnoise)
mixed = normalize(mixed)
plot_norm_compare_inverse_auto_correlation(mixed, normalize(gaussnoise), 16)
for ratio in np.linspace(0.0, 1.0, 10):
mixed = ratio * normalize(sindata) + (1.0 - ratio) * normalize(corr_gaussnoise)
mixed = normalize(mixed)
plot_norm_compare_inverse_auto_correlation(mixed, normalize(corr_gaussnoise), 16)
In [8]:
# 自己相関行列の逆の絶対値、対数絶対値の表示
def plot_abs_inverse_auto_correlation(data, dim):
S = calc_auto_corration_matrix(data, dim)
Sinv = np.linalg.inv(S)
Sinv /= np.linalg.norm(Sinv)
_, axs = plt.subplots(ncols=2, figsize=(10, 4))
sns.heatmap(np.abs(Sinv), ax=axs[0])
axs[0].set_title('Absolute Inverse Auto Correlation')
sns.heatmap(20 * np.log10(np.abs(Sinv)), ax=axs[1])
axs[1].set_title('Log-Absolute Inverse Auto Correlation [dB]')
for ratio in np.linspace(0.0, 1.0, 10):
mixed = ratio * normalize(sindata) + (1.0 - ratio) * normalize(gaussnoise)
mixed = normalize(mixed)
plot_abs_inverse_auto_correlation(mixed, 16)
for ratio in np.linspace(0.0, 1.0, 10):
mixed = ratio * normalize(sindata) + (1.0 - ratio) * normalize(corr_gaussnoise)
mixed = normalize(mixed)
plot_abs_inverse_auto_correlation(mixed, 16)
In [9]:
# Sherman–Morrison法による自己相関の逆の逐次計算
def calc_inv_auto_corr_by_sm(data, dim, forget):
Sinv = np.eye(dim)
for smpl in range(dim, data.size, 1):
xvec = data[smpl - dim:smpl].reshape((dim, 1))
Sx = np.dot(Sinv, xvec)
Sxxt = np.dot(Sx, Sx.T)
Sinv -= Sxxt / (forget + np.dot(xvec.T, Sx))
Sinv /= forget
return Sinv
# 自己相関行列の逆の絶対値、対数絶対値の表示
def plot_auto_corr_compare_sm_and_gt(data, dim):
S = calc_auto_corration_matrix(data, dim)
Sinv = np.linalg.inv(S)
Sinv_sm = calc_inv_auto_corr_by_sm(data, dim, 0.98)
Sinv /= np.linalg.norm(Sinv)
Sinv_sm /= np.linalg.norm(Sinv_sm)
_, axs = plt.subplots(ncols=3, figsize=(15, 4))
sns.heatmap(Sinv, ax=axs[0])
axs[0].set_title('Inverse Auto Correlation(Grand Trues)')
sns.heatmap(Sinv_sm, ax=axs[1])
axs[1].set_title('Inverse Auto Correlation(Sherman-Morrison)')
sns.heatmap(np.abs(Sinv - Sinv_sm), ax=axs[2])
axs[2].set_title('Absolute Error')
plot_auto_corr_compare_sm_and_gt(gaussnoise, 16)
plot_auto_corr_compare_sm_and_gt(normalize(corr_gaussnoise), 16)
# ノイズ以外の要素は、もはや気にしない!(上の考察より)
# plot_auto_corr_compare_sm_and_gt(normalize(sindata + gaussnoise), 16)
# plot_auto_corr_compare_sm_and_gt(normalize(sindata + corr_gaussnoise), 16)
In [10]:
from sklearn.covariance import GraphicalLasso
import scipy as sp
# Graphical Lassoによる精度行列計算(平均0だから自己相関行列の逆)
def calc_inv_auto_corr_by_gl(data, dim, alpha):
X = np.zeros((data.size - dim, dim))
for smpl in range(dim, data.size, 1):
X[smpl - dim:] = data[smpl - dim:smpl]
gl = GraphicalLasso(alpha=alpha)
gl.fit(X)
return gl.precision_
# 自己相関行列の逆の絶対値、対数絶対値の表示
def plot_auto_corr_compare_gl_and_gt(data, dim, alpha):
S = calc_auto_corration_matrix(data, dim)
Sinv = np.linalg.inv(S)
Sinv_gl = calc_inv_auto_corr_by_gl(data, dim, alpha)
Sinv /= np.linalg.norm(Sinv)
Sinv_gl /= np.linalg.norm(Sinv_gl)
_, axs = plt.subplots(ncols=3, figsize=(15, 4))
sns.heatmap(Sinv, ax=axs[0])
axs[0].set_title('Inverse Auto Correlation(Grand Trues)')
sns.heatmap(Sinv_gl, ax=axs[1])
axs[1].set_title('Inverse Auto Correlation(Graphical LASSO alpha: ' + str(alpha) + ')')
sns.heatmap(np.abs(Sinv - Sinv_gl), ax=axs[2])
axs[2].set_title('Absolute Error')
plot_auto_corr_compare_gl_and_gt(gaussnoise, 16, 0.1)
plot_auto_corr_compare_gl_and_gt(corr_gaussnoise, 16, 0.1)
考察¶
- 残差を見ると、Sherman-Morrisonよりも良さそうに見える。
- いよいよ、実装してみるべきなのかもしれない。
- 実データに対しても同様の傾向なのかは見ておきたい。
成果報告向けプロット¶
In [11]:
# 真のRの逆行列はどうなっているか?
plot_gt_auto_correlation(normalize(gaussnoise), 16)
plt.savefig("r_invr_gt_whitenoise.png")
plot_gt_auto_correlation(normalize(corr_gaussnoise), 16)
plt.savefig("r_invr_gt_corrnoise.png")
# Sherman-Morrisonとの比較
plot_auto_corr_compare_sm_and_gt(gaussnoise, 16)
plt.savefig("invr_comp_sm_and_gt_whitenoise.png")
plot_auto_corr_compare_sm_and_gt(normalize(corr_gaussnoise), 16)
plt.savefig("invr_comp_sm_and_gt_corrnoise.png")
# glassoとの比較
plot_auto_corr_compare_gl_and_gt(gaussnoise, 16, 0.1)
plt.savefig("invr_comp_gl_and_gt_whitenoise.png")
plot_auto_corr_compare_gl_and_gt(corr_gaussnoise, 16, 0.1)
plt.savefig("invr_comp_gl_and_gt_corrnoise.png")
追試:RLSの更新則で見ている自己相関行列¶
- Sharman-Morrisonで逐次計算しないやり方に該当する。観察。
- 性能としてはRLSの方が良いことが分かっている。
- 上とは何か違うものを見ている可能性が大。
In [12]:
# RLSが考えている自己相関行列の計算
def calc_inv_auto_corr_by_rls(data, dim, forget):
S = np.eye(dim)
for smpl in range(dim, data.size, 1):
xvec = data[smpl - dim:smpl].reshape((dim, 1))
S = forget * S + xvec @ xvec.T
# S = (smpl * S + xvec @ xvec.T) / (smpl + 1)
return S, np.linalg.inv(S)
# 自己相関行列の逆の絶対値、対数絶対値の表示
def plot_auto_corr_compare_rls_and_gt(data, dim, forget):
S = calc_auto_corration_matrix(data, dim)
Sinv = np.linalg.inv(S)
Srls, Sinv_rls = calc_inv_auto_corr_by_rls(data, dim, forget)
Sinv /= np.linalg.norm(Sinv)
Sinv_rls /= np.linalg.norm(Sinv_rls)
_, axs = plt.subplots(ncols=3, figsize=(15, 4))
# sns.heatmap(Sinv, ax=axs[0])
# axs[0].set_title('Inverse Auto Correlation(Grand Trues)')
sns.heatmap(Srls, ax=axs[0])
axs[0].set_title('Auto Correlation(RLS) Num Samples:%d' % data.size)
sns.heatmap(Sinv_rls, ax=axs[1])
axs[1].set_title('Inverse Auto Correlation(RLS)')
sns.heatmap(np.abs(Sinv - Sinv_rls), ax=axs[2])
axs[2].set_title('Absolute Error between the GT')
gaussnoise = np.random.normal(0, 1, NUM_SAMPLES)
corr_gaussnoise = make_1st_correlation(gaussnoise, 0.95)
for nsmpls in np.arange(100, 1000, 100):
plot_auto_corr_compare_rls_and_gt(gaussnoise[0:nsmpls], 16, 0.999)
plot_auto_corr_compare_rls_and_gt(corr_gaussnoise[0:nsmpls], 16, 0.999)
所感¶
- ほぼリファレンス(GT)と変わらないように見えるのだが...
- 各サンプルでGTに近ければ良い結果は出るはず。
- 定常状態(GTで得られる答え)との誤差は十分に小さいように見受けられる。
- 自分が考えてるAR(1)仮定をぶち込んで見てみよう。
In [13]:
# AR(1)仮定による逆行列計算
def calc_inv_auto_corr_by_ar1(data, dim, forget):
var00 = 1.0
var01 = 0.0
tri_mask = np.multiply(np.tri(dim, dim, -1), np.tri(dim, dim, 1).T)
tri_mask += tri_mask.T
# for smpl in range(dim, data.size, 1):
# var00 = forget * var00 + (1.0 - forget) * data[smpl] ** 2.0
# var01 = forget * var01 + (1.0 - forget) * data[smpl] * data[smpl - 1]
var00 = np.sum(data ** 2) / data.size
var01 = np.sum(data[0:data.size - 1] * data[1:data.size]) / (data.size - 1)
cor01 = var01 / var00
Rinv = (1.0 + (cor01 ** 2)) * np.eye(dim)
Rinv[0, 0] = Rinv[dim-1, dim-1] = 1.0
Rinv += -cor01 * tri_mask
Rinv /= var00 * (1.0 - (cor01 ** 2))
return Rinv
# 自己相関行列の逆の絶対値、対数絶対値の表示
def plot_auto_corr_compare_ar1_and_gt(data, dim, forget):
S = calc_auto_corration_matrix(data, dim)
Sinv = np.linalg.inv(S)
Sinv_ar1 = calc_inv_auto_corr_by_ar1(data, dim, forget)
Sinv_gl = calc_inv_auto_corr_by_gl(data, dim, 0.2)
Sinv /= np.linalg.norm(Sinv)
Sinv_ar1 /= np.linalg.norm(Sinv_ar1)
Sinv_gl /= np.linalg.norm(Sinv_gl)
_, axs = plt.subplots(ncols=3, figsize=(15, 4))
sns.heatmap(Sinv, ax=axs[0])
axs[0].set_title('Inverse Auto Correlation(Grand Trues)')
# sns.heatmap(Sinv_gl, ax=axs[0])
# axs[0].set_title('Inverse Auto Correlation(glasso)')
sns.heatmap(Sinv_ar1, ax=axs[1])
axs[1].set_title('Inverse Auto Correlation(AR(1))')
sns.heatmap(np.abs(Sinv - Sinv_ar1), ax=axs[2])
axs[2].set_title('Absolute Error')
for nsmpls in np.arange(30, 200, 20):
plot_auto_corr_compare_ar1_and_gt(gaussnoise[0:nsmpls], 16, 0.999)
plot_auto_corr_compare_ar1_and_gt(corr_gaussnoise[0:nsmpls], 16, 0.999)
In [29]:
# AR(2)仮定による逆行列計算
def calc_inv_auto_corr_by_ar2(data, dim, forget):
var00 = 1.0
var01 = 0.0
var02 = 0.0
tri_mask_ord1 = np.multiply(np.tri(dim, dim, -1), np.tri(dim, dim, 1).T)
tri_mask_ord1 += tri_mask_ord1.T
tri_mask_ord2 = np.multiply(np.tri(dim, dim, -2), np.tri(dim, dim, 2).T)
tri_mask_ord2 += tri_mask_ord2.T
# for smpl in range(dim, data.size, 1):
# var00 = forget * var00 + (1.0 - forget) * data[smpl] ** 2.0
# var01 = forget * var01 + (1.0 - forget) * data[smpl] * data[smpl - 1]
var00 = np.sum(data ** 2) / data.size
var01 = np.sum(data[0:data.size - 1] * data[1:data.size]) / (data.size - 1)
var02 = np.sum(data[0:data.size - 2] * data[2:data.size]) / (data.size - 2)
cor01 = var01 / var00
cor02 = var02 / var00
phi1 = cor01 * (1.0 - cor02) / (1.0 - (cor01 ** 2))
phi2 = (cor02 - (cor01 ** 2)) / (1.0 - (cor01 ** 2))
diag = (1.0 + (phi1 ** 2) + (phi2 ** 2)) * np.eye(dim)
diag[0, 0] = diag[dim-1, dim-1] = 1.0
diag[1, 1] = diag[dim-2, dim-2] = 1.0 + (phi1 ** 2)
tri_mask_ord1 = -phi1 * (1.0 - phi2) * tri_mask_ord1
tri_mask_ord1[0, 1] = tri_mask_ord1[1, 0] \
= tri_mask_ord1[dim-1, dim-2] = tri_mask_ord1[dim-2, dim-1] = -phi1
tri_mask_ord2 = -phi2 * tri_mask_ord2
Rinv = diag + tri_mask_ord1 + tri_mask_ord2
Rinv /= var00 * (1.0 + phi2) * (((1.0 - phi2) ** 2) - (phi1 ** 2)) / (1.0 - phi2)
return Rinv
# 自己相関行列の逆の絶対値、対数絶対値の表示
def plot_auto_corr_compare_ar2_and_gt(data, dim, forget):
S = calc_auto_corration_matrix(data, dim)
Sinv = np.linalg.inv(S)
Sinv_ar2 = calc_inv_auto_corr_by_ar2(data, dim, forget)
Sinv_gl = calc_inv_auto_corr_by_gl(data, dim, 0.1)
Sinv /= np.linalg.norm(Sinv)
Sinv_ar2 /= np.linalg.norm(Sinv_ar2)
Sinv_gl /= np.linalg.norm(Sinv_gl)
_, axs = plt.subplots(ncols=3, figsize=(15, 4))
sns.heatmap(Sinv, ax=axs[0])
axs[0].set_title('Inverse Auto Correlation(Grand Trues)')
# sns.heatmap(Sinv_gl, ax=axs[0])
# axs[0].set_title('Inverse Auto Correlation(glasso)')
sns.heatmap(Sinv_ar2, ax=axs[1])
axs[1].set_title('Inverse Auto Correlation(AR(2))')
sns.heatmap(np.abs(Sinv - Sinv_ar2), ax=axs[2])
axs[2].set_title('Absolute Error between the glasso')
# for nsmpls in np.arange(30, 200, 20):
for nsmpls in np.arange(100, 1000, 200):
plot_auto_corr_compare_ar2_and_gt(gaussnoise[0:nsmpls], 16, 0.999)
plot_auto_corr_compare_ar2_and_gt(corr_gaussnoise[0:nsmpls], 16, 0.999)
所感¶
- AR(1), AR(2)(今実装した)は誤差がそれなりに(0.02くらい)出ている。
- RLSの考え方、つまりSherman-Morrisonの誤差より大きい。
- スパースなglassoで求めた結果と比較しても傾向は同じ。
In [15]:
# 自己相関行列の逆の絶対値、対数絶対値の表示
def plot_auto_corr_compare_ar1_ar2_gt(data, dim, forget):
S = calc_auto_corration_matrix(data, dim)
Sinv = np.linalg.inv(S)
Sinv_ar1 = calc_inv_auto_corr_by_ar1(data, dim, forget)
Sinv_ar2 = calc_inv_auto_corr_by_ar2(data, dim, forget)
_, Sinv_rls = calc_inv_auto_corr_by_rls(data, dim, forget)
Sinv /= np.linalg.norm(Sinv)
Sinv_ar1 /= np.linalg.norm(Sinv_ar1)
Sinv_ar2 /= np.linalg.norm(Sinv_ar2)
Sinv_rls /= np.linalg.norm(Sinv_rls)
_, axs = plt.subplots(ncols=5, figsize=(16, 3))
sns.heatmap(Sinv_ar1, ax=axs[0])
axs[0].set_title('Inverse Auto Correlation(AR(1))')
sns.heatmap(Sinv_ar2, ax=axs[1])
axs[1].set_title('Inverse Auto Correlation(AR(2))')
sns.heatmap(np.abs(Sinv - Sinv_ar1), ax=axs[2], vmin=0, vmax=0.1)
axs[2].set_title('Absolute Error (AR(1))')
sns.heatmap(np.abs(Sinv - Sinv_ar2), ax=axs[3], vmin=0, vmax=0.1)
axs[3].set_title('Absolute Error (AR(2))')
sns.heatmap(np.abs(Sinv - Sinv_rls), ax=axs[4], vmin=0, vmax=0.1)
axs[4].set_title('Absolute Error (RLS)')
In [16]:
# i.i.d. に対して実験
for nsmpls in np.arange(50, 600, 50):
print(nsmpls)
plot_auto_corr_compare_ar1_ar2_gt(gaussnoise[0:nsmpls], 16, 0.999)
plt.show()
In [17]:
# 相関付きに対して実験
for nsmpls in np.arange(50, 600, 50):
print(nsmpls)
plot_auto_corr_compare_ar1_ar2_gt(corr_gaussnoise[0:nsmpls], 16, 0.999)
plt.show()
In [18]:
# 真の逆行列との誤差RMSEの変化をプロット
def plot_inv_auto_corr_rmse_curve(dim, forget, smpl_range, data_generator, num_trials):
num_ranges = smpl_range.size
list_ar1 = np.zeros(num_ranges)
list_ar2 = np.zeros(num_ranges)
list_rls = np.zeros(num_ranges)
list_gl = np.zeros(num_ranges)
for _ in range(num_trials):
for idx, nsmpl in enumerate(smpl_range):
data = data_generator(nsmpl)
S = calc_auto_corration_matrix(data, dim)
Sinv = np.linalg.inv(S)
Sinv_ar1 = calc_inv_auto_corr_by_ar1(data, dim, forget)
Sinv_ar2 = calc_inv_auto_corr_by_ar2(data, dim, forget)
_, Sinv_rls = calc_inv_auto_corr_by_rls(data, dim, forget)
Sinv_gl = calc_inv_auto_corr_by_gl(data, dim, 1e-1)
Sinv /= np.linalg.norm(Sinv)
Sinv_ar1 /= np.linalg.norm(Sinv_ar1)
Sinv_ar2 /= np.linalg.norm(Sinv_ar2)
Sinv_rls /= np.linalg.norm(Sinv_rls)
Sinv_gl /= np.linalg.norm(Sinv_gl)
list_ar1[idx] += np.linalg.norm(Sinv - Sinv_ar1)
list_ar2[idx] += np.linalg.norm(Sinv - Sinv_ar2)
list_rls[idx] += np.linalg.norm(Sinv - Sinv_rls)
list_gl[idx] += np.linalg.norm(Sinv - Sinv_gl)
list_ar1 /= num_trials
list_ar2 /= num_trials
list_rls /= num_trials
list_gl /= num_trials
plt.plot(smpl_range, list_ar1, label='AR(1)')
plt.plot(smpl_range, list_ar2, label='AR(2)')
plt.plot(smpl_range, list_rls, label='RLS(S-M)')
plt.plot(smpl_range, list_gl, label='glasso')
plt.legend()
plt.xticks(smpl_range[::num_ranges//10])
# plt.yscale('log')
plt.grid()
plt.show()
def generate_normal(num_samples):
return np.random.normal(0, 1, num_samples)
def generate_corr_normal(num_samples):
return normalize(make_1st_correlation(generate_normal(num_samples), 0.95))
In [19]:
# i.i.d.雑音に対して試す
plot_inv_auto_corr_rmse_curve(4, 0.999, np.arange(50, 1000, 20), generate_normal, 20)
plot_inv_auto_corr_rmse_curve(16, 0.999, np.arange(50, 1000, 20), generate_normal, 10)
In [20]:
# 相関付き雑音に対して試す
plot_inv_auto_corr_rmse_curve(4, 0.999, np.arange(50, 1000, 20), generate_corr_normal, 20)
plot_inv_auto_corr_rmse_curve(16, 0.999, np.arange(50, 1000, 20), generate_corr_normal, 10)
考察¶
- AR(1), AR(2)は大数の法則的に漸近しているように見える(1/√n)。
- これは最尤推定量を使ったからと想像。指数移動平均でやっても結果はほぼ同じだった。
- これだと遅い、ということなのかもしれない。
- 最初の推定値が悪すぎるらしい。それによってRLSと差が出てしまっているっぽい。
- RLSの推定値が "良い" とはどういうことなのだろうか。。。
- AR(2)の方がわずかに残差が小さい傾向。
- glassoはARよりちょっと良いくらいか?相関付き雑音だと却って悪いくらいの性能。
- 正則化パラメータを小さくすると良くなるが、相関付き雑音で条件数が悪くなって警告・エラーを出し始める
次確かめること¶
- 大域的な真値への収束特性確認。AR(1)が良かった記憶があるが、果たして本当にそうだったか。
In [21]:
# 真の逆行列(固定)との誤差RMSEの変化をプロット
def plot_fixed_inv_auto_corr_rmse_curve(dim, forget, smpl_range, data_generator, num_trials):
num_ranges = smpl_range.size
list_ar1 = np.zeros(num_ranges)
list_ar2 = np.zeros(num_ranges)
list_rls = np.zeros(num_ranges)
list_gl = np.zeros(num_ranges)
for _ in range(num_trials):
# 長期サンプルにおける真の逆行列を先に求めておく
wholedata = data_generator(np.max(smpl_range))
S = calc_auto_corration_matrix(wholedata, dim)
Sinv = np.linalg.inv(S)
Sinv /= np.linalg.norm(Sinv)
for idx, nsmpl in enumerate(smpl_range):
data = wholedata[0:nsmpl]
Sinv_ar1 = calc_inv_auto_corr_by_ar1(data, dim, forget)
Sinv_ar2 = calc_inv_auto_corr_by_ar2(data, dim, forget)
_, Sinv_rls = calc_inv_auto_corr_by_rls(data, dim, forget)
# Sinv_gl = calc_inv_auto_corr_by_gl(data, dim, 1e-1)
Sinv_ar1 /= np.linalg.norm(Sinv_ar1)
Sinv_ar2 /= np.linalg.norm(Sinv_ar2)
Sinv_rls /= np.linalg.norm(Sinv_rls)
# Sinv_gl /= np.linalg.norm(Sinv_gl)
list_ar1[idx] += np.linalg.norm(Sinv - Sinv_ar1)
list_ar2[idx] += np.linalg.norm(Sinv - Sinv_ar2)
list_rls[idx] += np.linalg.norm(Sinv - Sinv_rls)
# list_gl[idx] += np.linalg.norm(Sinv - Sinv_gl)
list_ar1 /= num_trials
list_ar2 /= num_trials
list_rls /= num_trials
# list_gl /= num_trials
plt.plot(smpl_range, list_ar1, label='AR(1)')
plt.plot(smpl_range, list_ar2, label='AR(2)')
plt.plot(smpl_range, list_rls, label='RLS(S-M)')
# plt.plot(smpl_range, list_gl, label='glasso')
plt.legend()
plt.xticks(smpl_range[::num_ranges//10])
# plt.yscale('log')
plt.grid()
plt.show()
In [22]:
# i.i.d.雑音に対して試す
plot_fixed_inv_auto_corr_rmse_curve(4, 0.999, np.arange(20, 700, 10), generate_normal, 20)
plot_fixed_inv_auto_corr_rmse_curve(8, 0.999, np.arange(20, 700, 10), generate_normal, 20)
plot_fixed_inv_auto_corr_rmse_curve(16, 0.999, np.arange(20, 700, 10), generate_normal, 20)
In [23]:
# 相関付き雑音に対して試す
plot_fixed_inv_auto_corr_rmse_curve(4, 0.999, np.arange(20, 700, 10), generate_corr_normal, 20)
plot_fixed_inv_auto_corr_rmse_curve(8, 0.999, np.arange(20, 700, 10), generate_corr_normal, 20)
plot_fixed_inv_auto_corr_rmse_curve(16, 0.999, np.arange(20, 700, 10), generate_corr_normal, 20)
報告用プロット¶
In [24]:
plot_auto_corr_compare_ar1_and_gt(gaussnoise, 16, 0.999)
plt.savefig("InverseAutoCorrCompareAR1_iid.png")
plot_auto_corr_compare_ar1_and_gt(corr_gaussnoise, 16, 0.999)
plt.savefig("InverseAutoCorrCompareAR1_corr.png")
追加実験: 実データの自己相関行列の逆¶
In [54]:
import soundfile
realdata, _ = soundfile.read('./data/seijya_no_kousin_head4s_Lch.wav')
for smpl in np.arange(5000, 50000, 5000):
partial_data = realdata[smpl-5000:smpl]
print('%d-%d samples:' % (smpl-5000, smpl))
plot_auto_corr_compare_ar1_and_gt(partial_data, 16, 0.999)
plot_auto_corr_compare_ar2_and_gt(partial_data, 16, 0.999)
plot_auto_corr_compare_gl_and_gt(partial_data, 16, 0.01)
plt.show()
考察¶
- 実データの自己相関行列がほぼ変わらない
- 真値:行列の中央が最大値(約0.2)で、端に行くほど減衰。
- 真値:微妙に波打ってる。
- 真値、AR(1)はずっと変わってない
- (1,1), (N,N)で大きな差異が出ているのが気になる。
- glassoは僅かに変わっているが、対角優位で端に行くほど減衰する傾向は変わらない
- AR(1):glassoの推定結果に近い
In [ ]: